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ABSTRACT 

Context. Observations have long demonstrated the molecular diversity of the diffuse interstellar medium (ISM). Only now, with the 
advent of high-performance computing, does it become possible for numerical simulations of astrophysical fluids to include a treat- 
ment of chemistry, in order to faithfully reproduce the abundances of the many observed species, and especially that of CO, which is 
used as a proxy for molecular hydrogen. When applying photon-dominated region (PDR) codes to describe the UV-driven chemistry 
of uniform density cloud models, it is found that the observed abundances of CO are not well reproduced. 

Aims. Our main purpose is to estimate the effect of assuming uniform density on the line-of-sight in PDR chemistry models, compared 
to a more realistic distribution for which total gas densities may well vary by several orders of magnitude. A secondary goal of this 
paper is to estimate the amount of molecular hydrogen which is not properly traced by the CO (/ = 1 — > 0) line, the so-called "dark 
molecular gas". 

Methods. We use results from a magnetohydrodynamical (MHD) simulation as a model for the density structures found in a turbulent 
diffuse ISM with no star-formation activity. The Meudon PDR code is then applied to a number of lines of sight through this model, 
to derive their chemical structures. 

Results. It is found that, compared to the uniform density assumption, maximal chemical abundances for H 2 , CO, CH and CN 
are increased by a factor ~ 2 - 4 when taking into account density fluctuations on the line of sight. The correlations between 
column densities of CO, CH and CN with respect to those of H 2 are also found to be in better overall agreement with observa- 
tions. For instance, at N(H 2 ) > 2 10 20 cm" 2 , while observations suggest that d[log Af(CO)]/d[log N(H 2 )] - 3.07 + 0.73, we find 
d[logA^(CO)]/d[logA^(H 2 )] - 14 when assuming uniform density, and d[log iV(CO)]/d[log N(H 2 )] - 5.2 when including density 
fluctuations. 

Key words. ISM : structure - ISM : clouds - ISM : photon-dominated region (PDR) 



1. Introduction 

The interstellar medium (ISM) is a complex system : its struc- 
ture and dynamics are governed by the interaction of many pro- 
cesses covering a wide range of scales, involving both micro- 
and macrophysics. To understand how the ISM works is also 
essential on the pathway to star and planetary system forma- 
tion. This field has seen much progress in recent years, both on 
the observational side, with the results from the Herschel Space 
— 20l0||de Graauw et al. ||20l0^ , and 



Observatory ( Pilbratt et al. 



on the theoretical side, with ever-improving numerical simula- 
tions (see e.g. Banerjee et al. | 2009 1 |Vazquez-S emadeni et al. 



2010), which now consistently treat self-gravity, thermodynam- 
ics and magnetohydrodynamics (MHD). 

The new challenge for numerical simulations of the ISM 
is now to incorporate some treatment of chemistry and grain 
physics, in order to compare with observations of atomic and 
molecular lines. One possibility, followed for instance by |Glover 
et al. | ([2010) is to perform multi-fluid simulations, each fluid 



representing a chemical species and being coupled to the others 
via a network of reactions. This approach has the advantage that 
it naturally solves for the time-dependent distribution and dy- 
namics of the various species in three-dimensional space, but its 
computational cost requires making use of simplifying assump- 
tions, mo st notably a sm all number of species and reactions. For 
instance, [Glover et al. | ( |2010| ) treat a network of 218 reactions 



for 32 species, 13 of which are assumed to be in instantaneous 
chemical equilibrium, so that only 19 are actually fully treated 
as time-dependent quantities. 

In this paper, we follow a different approach : we post- 
process a single-fluid MHD simulation with the Meudon PDR 
code. This has the advantage of providing a full chemical net- 
work (99 species, 1362 reactions), but on the other hand, it im- 
plies a one-dimensional, steady-state treatment. With this ap- 
proach, our main goal is to discuss the effects of realistic den- 
sity fluctuations on PDR chemistry, since previous studies have 



focused on uniform density models, as in |Le Petit et al. |([20 06), 
or on simplified clumpy models, such as those by Wolfire et al. 



( |2010| ). A second goal of this study is to estimate the amount of 
"dark molecular gas" that can be expected from observations of 
the diffuse ISM, i.e. gas where hydrogen is mostly in the form of 
H2, but where CO is too scarce to be seen in the usual tracer that 
is the (/ = 1 — > 0) line, given current sensitivities. We wish to 



compare this estimate with available observations ( Grenier et al. 


2005 , Leroy et al. 2007 Abdo et al. 


2010( |Velusamy et al. 


2010) and models flWolfire et al. ||2010 





The paper is organised as follows : Section [2] describes the 
MHD simulation used and presents the PDR code in its current 
state, while section [3] gives an overview of the post-processing 
method used to combine the two. Section [4] discusses the lines 
of sight selected for running our analysis. Results and com- 
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parison to observations are presented in section [5] Section [6] 
offers further discussion and a summary. Details on using the 
Meudon PDR code with density profiles and our strategy for 
post-processing results can be found in the appendices at the end 
of the paper. 

2. Tools 

2.1. The MHD simulation 

As a model of the turbulent diffuse interstellar medium, we use 
data from a magnetohydrodynamical (MHD) simulation per- 
formed by IHennebelle et al. | (|2008|) u sing the RAMSES code 
(Teyssier , 2002; Fromang et al. 2006 ). A notable advantage 
of this code resides in its adaptive mesh refinement (AMR) ca- 
pabilities, making it able to locally reach extremely high spa- 
tial resolutions. The initial setup of the simulation is a cube of 
homogeneous warm neutral gas (WNM) of atomic gas (a 10:1 
mixture of H and He) with density nn = 1 cm -3 and temper- 
ature T = 8000 K. The cube is L = 50 pc on each side, and 
two converging flows of WNM are injected from opposing faces 
along the X axis with relative velocity AV X ~ 40 km.s -1 , which 
means that the Mach number of each flow with respect to the 
ambient WNM is Al ~ 2. Transverse and longitudinal velocity 
modulations are imposed, with amplitudes roughly equal to the 
mean flow velocity (20 km.s -1 ). Periodic boundary conditions 
are applied on the remaining four faces. The simulation starts on 
a 25 6 3 grid, with two extra levels of refinement based on den- 
sity thresholds, so that the effective resolution of the simulation 
is ~ 0.05 pc. For our purposes, we regridded the data regularly 
on a 1024 3 cube, therefore using the maximum resolution over 
the whole domain. A magnetic field is present in the simula- 
tion, and is initially parallel to the X axis with an intensity of 
abo ut 5 //G, consistent with observational values at these densi- 
ties (Crutcher et al. | [2010). The converging flows collide near 
the midplane X = about 1 Myr into the simulation, and as the 
total mass grows from ~ 3000 M initially to about 10 times 
this value at the end of the simulation (t = 13.11 Myr), gravity 
eventually takes over and, combined with the effects of thermal 
instability ([Field] [T9651 [Hennebelle & Perault ||2 000; Koyama| 
& Inutsuka , 2002; Heitsch ^eTaT. ||2005t IHennebelle & Audit [ 



2007 ; Banerj ee et al. ||2009| ), leads to the formation of cold dense 
clumps (nn > 100 cm~\ while T ~ 10-50 K) within a much 
more diffuse and warm interclump medium (nn ~ 1 — 10 cm" 3 
and T ~ 10 3 - 10 4 K). This occurs at t ~ 12 Myr. See|Hennebelle 



et al. ( 2008) for a more detailed description. Fig. [1] shows the 



column density of the gas viewed along the X axis, which is that 
of the incoming flows. 

2.2. The Meudon PDR code 



The Meudon PDR code (ht tp://pdr.obspm.fr/) is a pub- 
licly avail able set of routines (|Le Bourlot et aI~[|1993[|Le Petit et| 
200 6; Gonzalez Garci a et al. | 2008| ) whose purpose is to de- 



al. 



scribe the UV-driven chemistry of interstellar clouds, and in par- 
ticular of photon-dominated regions (PDRs). It is a steady- state 
one-dimensional code in which a plane-parallel slab of gas and 
dust is illuminated on either or both sides by the light from a star 
or by the standard Inter Stellar Radiation Field (ISRF), which is 
d efined using expressions from Mat his et al. | ( |1983] ) and [Black 
[ (1994). Usually run with homogeneous gas densities, the code 
can also accept density profiles input via a text file. At each point 
along the line of sight, radiative transfer in the UV is treated to 
solve for the H/H2 transition, using either the approximations by 



|Federman et al. | ( [1979| ), or an exact method based o n a spheri- 
cal harmonics expansion of the specific intensity ( Goicoechea & 
|Le Bourlot] [2007). A number of heating (photoelectric effect on 
grains, cosmic rays) and cooling (infrared and millimeter emis- 
sion lines) processes contribute to the computation of thermal 
balance. Outputs of the code include gas properties such as tem- 
perature and ionization fraction, radiation energy density, chem- 
ical abundances and column densities, level populations and line 
intensities. The code is iterative and therefore requires the user 
to check whether convergence has been achieved. 

3. Method overview 

Naively, one would like to run the PDR code on all lines of sight 
through the simulated cube to derive a three-dimensional chem- 
ical structure, as well as line-of-sight integrated observables (an 
emission map for the CII [158 fim] line, for instance). However, 
this is neither feasible nor desirable. 

Two computational reasons preclude this brute force ap- 
proach. Firstly, the cost is prohibitive : just one global iteration 
of the PDR code for a single run typically completes in a cou- 
ple hours on a GNU/Linux machine equipped with 4 dual-core 
64-bit x86 processors and 64GB of memory. As a run usually 
converges in some 10 iterations, it takes about a day to process 
a single line of sight. The PDR code, however, has been ported 
on the EGEE gricQ which allows us to perform ~ 100 runs si- 
multaneously in the same timeframe. Still, it is not reasonable to 
consider treating more than ~ 10 3 lines of sight in this work. 

The second computational issue to consider is that the PDR 
code has trouble converging in low-density regions, of which 
there are many in the MHD simulation. Such are the vicinities 
of X = ±25 pc, where WNM gas (n H = 1 cm -3 ) is entering 
the box, but low density regions are actually found everywhere 
throughout the cube (the volume filling factor of regions where 
n H < 20 cm -3 is f v - 0.96). This convergence issue might be 
related to the fact that Ly-a emission, which is a major cooling 
process for nu < 5 cm -3 , is not included in the code, making 
the computation of thermal balance by the PDR code in these 
regions unreliable. 

Besides these computational hurdles, it should be noted that 
the code is one-dimensional, and treats a density profile as that 
of a plane-parallel slab of gas. Imagine then that a line of sight 
intercepts a dense clump of matter : the gas lying directly be- 
hind that clump will be shielded from incoming radiation, since 
the most energetic UV photons will have been absorbed by the 
clump. This leads to shadowing artifacts in regions which in re- 
ality may well be illuminated from other directions, as the ISM 
has a complex, fractal-like, and evolving structure. 

In this paper, we deal with these artifacts in the following 
way : the physical conditions and chemical composition at every 
grid point (X, 7, Z) considered in the analysis may be obtained 
by running the PDR code in p directions going through that 
grid point, for instance along the main coordinate axes. Thus, 
each quantity F(X, Y, Z) output by the code has p possible val- 
ues fi,fz, . . . ,f p . This is the case, in particular, of the radia- 
tion energy density E, which has possible values e\, £2, • • • , e p . 
To combine the p runs at this grid point, we select direction 
Po for which the radiation energy density there is maximum, 
e Po = max{ei}i <i<p . This choice is discussed in section [6] Of 
course po is a function of (X, 7, Z). Once this is done, we choose 
F(X, Y, Z) = f po for every quantity F output by the code. This 
procedure takes better account of the porosity of the simulated 
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Fig. 1. Total gas column density along the X axis of the MHD 
simulation snapshot used here. The box in the upper left side 
marks the position of the clump selected for running our analy- 
sis. 



structures to the ISRF, while ensuring element conservation at 
each grid point. Ideally, the more directions p, the better, but 
this obviously comes with an increased computational cost. In 
this paper, we choose p = 2 as a compromise, running the PDR 
code in two orthogonal directions. 

The analysis in this paper focuses on a small subset of a sin- 
gle simulation snapshot, which we discuss in the next section. 
After applying a density threshold no = 20 cm -3 , which is the 



lowest value considered in the grid of models run by Le Petit 
|et al. | ([2006 ) and therefore deemed sufficient to ensure conver- 
gence, we extract a number of one-dimensional density profiles 
from this thresholded subset. We run the PDR code on these pro- 
files and combine results to derive the chemical structure of the 
clump. 



5.0-24.5-24.0-23.5-23.0-22.5-22.0-21.5-21.0 

npc] 



Fig. 2. Maximum gas density nn along the lines of sight paral- 
lel to the X axis, within the selected clump. Contours show the 
total gas column density Nh from 3 10 21 cm -2 to 1.1 10 22 cm -2 
in steps of 10 21 cm -2 . The "clump" displayed here has a size 
roughly 1.5 pc x 2.5 pc. The 2D slab of gas under study is seen 
projected as the single-pixel- wide AB strip. The white triangle 
marks the position of the maximum value of max (nn) in this re- 
gion, and the white circle that of the maximum of Nh- They are 
separated by 0.65 pc. 



Table 1. Properties of the selected observational clump. (F) 
refers to the direct average and F to the density-weighted av- 
erage of any quantity F in this table. The mass is computed as 
M - "V/MiH (n H ), where *V is the 3D volume corresponding to 
the 2D clump, m# = 1.66 10" 24 g is the mass of the hydrogen 
atom and \i-\A corresponds to a 1 : 10 number ratio for He with 
respect to H. The turbulent velocity dispersion includes all three 
velocity components, <x 2 D = <x 2 + <x 2 + cr|. 



4. Lines of sight selected for PDR computations 

4.1. Simulation snapshot 

The snapshot chosen to run our analysis on is timed at 7.35 Myr, 
when the densest parts of the cloud reach ft max ~ 9.10 3 cm -3 . 
Some of the structures present in the simulation at this time are 
self-gravitating, but we are confident that they are still diffuse 
enough that the simulation snapshot is representative of a non- 
starforming region of the ISM. We may therefore run our analy- 
sis in the absence of any illuminating star, with the ISRF being 
the only source of primary UV photons. 



4.2. Selected clump 

To select a representative subset, we may not e that observationa l 
PDRs such as the Horsehead Nebula (see e.g. Pety et al. (2007)) 
are found at the edge of dense and cold clouds of gas and dust, 
illuminated by ambient FUV light and possibly nearby young 
stars. It thus makes sense to focus on a "clump", defined obser- 
vationally as a connected structure with a significantly higher 
column density than its surroundings. We identify clumps via a 
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friend-of-friend algorithm on the column density map along the 
X axis (Fig.[T]), using a threshold No = 3.10 21 cm -2 , which cor- 
responds to a mean density (nn) = no = 20 cm -3 over 50 pc. 

The selected clump, which lies at the top left corner of 
the simulation's field, harbours an interesting feature, shown on 
Fig. [2] That figure represents, in colour scale, the map of the 
maximum gas density max (nn) encountered along the X axis, 
for every line of sight within the clump. It so happens that the 
peak of N H does not match that of max (n H ), or even a local max- 
imum of the latter. This is important to note for species which 
may be sensitive to the local gas density rather than to the total 
column density. Properties of that selected clump are listed in 
Tableffl 
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Table 2. Properties of the 156 one-dimensional profiles extracted 
parallel to the X axis. Listed are the minimum, maximum and en- 
semble average values for the size, average density, column den- 
sity, visual extinction (corresponding to Nh via Eq.[T]), density- 
weighted average temperature and line-of-sight velocity disper- 
sion. For that last quantity, the minimum value is not shown, as 
it is too small to be meaningful. 



Parameter (F) 


min (F) 


max (F) 


(F) 


Size L [pc] 


0.15 


11.2 


2.18 


Density (n H ) [cm -3 ] 


20 


571 


155 


Column density N H [10 20 cm -2 ] 


0.117 


107 


11.9 


Visual extinction A v 


6.3 10- 3 


5.7 


0.64 


Temperature T [K] 


22 


924 


88 


Velocity dispersion <r x [km.s -1 ] 




1.9 


0.50 



4.3. Selected lines of sight 

The observational clump shown on Fig.[2]has a size AY x AZ ^ 
1.5 pc x 2.5 pc. In the X direction, most of its gas is located 
in the central AX ^15 pc around X = 0. Given the pixel size 
6 ^ 0.05 pc, applying the PDR code on a structure of that size 
along the three coordinate axes requires some 25000 runs, which 
is beyond the scope of this work. Consequently, we restrict our 
study to a 2D slab of gas across the observational clump. It is 
projected on Fig. [2] as the one-pixel-wide strip AB, which is 
therefore ~ 3.7 pc long, ~ 0.05 pc wide and contains 76 lines 
of sight along the X axis. These sample a wide range of col- 
umn densities, from Nh = 6.11 10 20 cm -2 (near the A end) 
to Nh - 1.11 10 22 cm -2 , corresponding to visual extinctions 
Ay = 0.33 to Ay = 5.93, using the conversion from total hydro- 
gen column density N H = N(H) + 2N(H 2 ) 

Rv(_N H 

Cr 



A v = 



,d V 1 cm 1 1 



(1) 



with 7? v = 3.1 andC D = 5.8 x 10 21 cnT^mag -1 (see Table 



A.l 



and |Le Petit et al. |[2006| ). The structure of the gas on these lines 
of sight (Fig. [3]) is complex, with many small, dense and cold 
regions {n H - 10 3 cm -3 , T - 20 K) interconnected via a fil- 
amentary structure and embedded within a much more diffuse 
and warm medium (n H - 1 cm -3 , T - 10 4 K). The overdense 
regions (uh > no) are located near the midplane of the simula- 
tion (-9 pc < X < 2 pc), where the flows collide and the gas 
condenses into cold structures. 

It is this subset of the simulated cube, shown on Fig. [3] 
which is the focus of our study, and from which we extract one- 
dimensional density profiles. As Fig.|4]shows, this subset indeed 
contains most of the gas on the lines of sight within the AB strip : 
Except in the outermost regions where Nh < 1.5 10 21 cm -2 , 
column densities along X over that region represent more than 
half the total column densities over the full 50 pc lines of sight. 
Fig.|4]also shows that, whithin the AB strip, the peaks of Nh and 
max (n H ) are still separated, by about 0.55 pc. 

As explained in section [3] we consider two (p = 2) possi- 
ble directions for the one-dimensional density profiles extracted 
from the 2D subset, namely those parallel to the X or Y axis. The 
dashed line on Fig.[3]marks the location of such an extracted pro- 
file, which is shown on Fig. [5] Note that we consider profiles to 
be connectedly overdense, which means that, along each line of 
sight, we may extract several profiles separated by underdense 
regions nn < no. Such is the line of sight parallel to the X axis 
located at Y = -23 pc, for instance. 

We thus extract 447 density profiles, of which 156 are paral- 
lel to the X axis and 291 are parallel to the Y axis. Their statis- 
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Fig. 3. Structure of the gas in the 2D slab under study. This is 
a close-up view on the region X ^ where the incoming flows 
of WNM collide to form cold structures, and only regions where 
rin > no = 20 cm -3 are shown. The colour image shows the 
gas density nn in logarithmic scale, while contours show the gas 
temperature T M hd at 20 K (higher densities), 30 K and 50 K 
(lower densities). The 2D projection of the velocity field (Vx, Vy) 
is also shown as yellow arrows, the length of which indicate the 
velocity modulus at that location (at the center of each arrow). 
The A and B extremities of the observational ID strip are indi- 
cated for reference, and the dashed line marks an example loca- 
tion for the extracted density profiles on which the PDR code is 
run (Fig. [5]). The grey areas are outside of the domain used for 
PDR computations. 
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Fig. 4. Total gas column densities along the X axis within the AB 
strip under study (top plot) and maximum gas density max (jih) 
on the same lines of sight (bottom plot). Shown are the col- 
umn densities over the full 50 pc lines of sight along the X axis 
(dashed line) and the column densities for the overdense regions 
n H > no shown on Fig. [3] (solid line). Grey areas mark lines of 
sight for which less than 50% of the mass is in the overdense 
region. The dash-dotted lines mark the positions of the maxima 
of N H and max (n H ), which are separated by 0.55 pc. 




Fig. 5. Example of a density profile used in the PDR code. This 
is the profile extracted at the location of the dashed line on Fig. [3] 

Table 3. Same as Table § but for the 291 one-dimensional pro- 
files extracted parallel to the Y axis. 



Parameter (F) 


min (F) 


max (F) 


(F) 


Size L [pc] 


0.24 


2.78 


1.17 


Density (n H ) [cm -3 ] 


22 


655 


188 


Column density N H [10 20 cm" 2 ] 


0.165 


31.6 


6.39 


Visual extinction A v 


8.8 10- 3 


1.7 


0.34 


Temperature T [K] 


21 


464 


56 


Velocity dispersion <r Y [km.s -1 ] 




1.51 


0.37 



Table 4. Typical threshold abundances X a used in Figs. 0and[8l 



C + : 10" 5 
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tical properties are summarized in Tables [2] and [3] emphasizing 
the large dynamic range they sample in column densities (~ 10 4 ) 
and mean densities (~ 30). The large values in density-weighted 
temperatures correspond to density profiles that never much de- 
viate from no = 20 cm -3 . For each of these 447 profiles, we 
run the PDR code assuming identical illumination on both sides. 
Since one of the objectives of this paper is to assess the effect 
of realistic density distributions along the line of sight on the 
chemical composition of interstellar clouds, we also apply the 
PDR code on a homogeneous reference model for each extracted 
profile. We specify this model, called uniform in the following, 
as having the same mean density (hh) and total visual extinc- 
tion Ay as the inhomogeneous model, which we dub los from 
now on. Unless otherwise specified, results presented in the next 
section refer to these los models. The setup for all runs is de- 
tailed in appendix [A] and their post-processing is described in 
appendix [B] 

5. Results 

5.1. Temperature comparison 

The PDR code and the MHD simulation both treat thermal bal- 
ance, so that we have two estimates of the gas temperature, 
which we can compare : Fig. [6] shows the ratio r = r PDR /r M HD 
of the gas temperature 7>dr output by the PDR code to the 
gas temperature Tmhd computed in the MHD simulation, at ev- 
ery point in the subset under study, versus the total gas density 
nu at that point. Average ratios (r) in selected density bins are 
also shown. It appears quite clearly that (r) is close to 1, with 
0.3<(r)<2.0 over the whole range of densities. 

The fact that Tpdr ~ ?mhd actually comes as a pleasant sur- 
prise, considering the differences between the PDR and MHD 
computations : while the former is ID, steady-state, and includes 
cooling via the infrared and submillimeter lines from atomic 
and molecular species, especially H 2 transitions ( |Le Petit et al. 
[| |2006| ), the latter is 3D, dynamical, and only includes cooling 
via the fine structure lines of CII [158/mi] and OI [63//m], as 
well as the recombination of electrons with ionized PAHs. This 
suggests that unless a very precise knowledge of the tempera- 
ture is needed, it is probably not necessary, at least as a first 
approximation and in the range of densities and temperatures 
probed here, to refine the details of cooling processes in our 
MHD simulations, as the simple cooling function currently used 
already yields gas temperatures close to those found using the 
more detailed p rocesses of the PDR code. Similar conclusions 
were reached by Glover et al. ( 2010| ). 



5.2. Chemical structure and comparison to observations 

The spatial distributions of H, H 2 , C + , C, CO, CS, CH and CN 
in the simulation subset are shown on Figs. [7] and [8] In these 
figures, points were clipped where the density n(a) of species a 
was below no, a = X a no, with X a a typical threshold abundance 
for the detection of that species (see Table [4]). Although shadow- 
ing effects remain (for instance on the atomic hydrogen map), 
these figures show how some species (e.g. CO, CN) trace denser 
gas than others (C + , CH). This appears more clearly when plot- 
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Fig. 6. Ratio of the gas temperatures Tpdr / ^mhd versus total 
gas density nn- Grey crosses show all points, while black cir- 
cles are average values over density bins in logarithmic scale, 
with error bars standing for +lcr. The dashed line corresponds 
to TVdr/^mhd = I- 



ting these abundances, averaged over density bins, versus total 
gas density hh (Fig. [5]) : C + traces gas uniformly up to hr > 
10 3 cm -3 , while CO starts rising up at n H > 250 cm -3 , right 
about where CH flattens out. This break in the slope for CO oc- 
curs after the molecular transition (fn 2 ) = 2(n(H2)/w#) = 1/2, 
which is at n H 100 cm -3 . Considering the abundances of C 
and C + , this means that a significant fraction of the molecular 
gas (i.e. where hydrogen is mostly in the form of H2) is bet- 
ter traced by C and C + than by CO. This "dark molecular gas" 
fraction is the sub ject of subsection |5.4| CH, on the other hand, 
nicely follows H2 ( Sheffer et al. 2008 ). To complete the picture, 
C and CN have a similar slope throughout the density range, al- 
though CN seems to break away slightly at n H > 10 3 cm -3 , to 
follow CO. 

To be more precise on these apparent correlations (CH vs. 
H 2 , CS vs. C and CN vs. CO), we show, on Fig. [10| abun- 
dance ratios for these similarly distributed species in the region 
where they are all significantly present, that is the cloudlet at 
(X -4.7 pc, Y - -24 pc) (see Figs. [7] and [8]). We can see that 
n(CU)/n(U 2 ) and n(CN)/n(CO) have a similar "ringlike" be- 
haviour, rising to maximum values ft(CH)/?z(H 2 ) ~ 5.5 10" 8 and 
n(CN)/n(CO) ~ 10" 2 at total gas densities n H ~ 400-500 cm" 3 , 
then falling at higher densities, to fz(CH)/^(H 2 ) ~ 4 10" 8 and 
7i(CN)/n(CO) ~ 2 10" 3 , respectively. For n(CS)/n(C), there is 
a slight loss of azimuthal symmetry, but the overall trend is the 
same, although maximum values of ~ 10" 3 are reached at larger 
densities nn ~ 1000 cm -3 before falling to ~ 3 10" 4 at the peak. 

The comparison to observational data requires us to work 
not with densities but with column densities, which are what ob- 
servers have access to. To his end, we compute column densities 
for CO, CH and CN on each line of sight parallel to the X or 
Y axis and plot them versus those of H 2 ( Fig. 
observational data plots in Sheff er et al. | (f2008 ) 



1 1 ), to match the 
. We use a colour 

scheme to specify the mean gas density (nn) on the line of sight, 
and plot separately data points corresponding to the long lines of 
sight parallel to X (squares) and to the shorter lines of si ght par- 
allel to Y (circles). Fits to observational data derived by Sheffer 
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Fig. 7. H (top left), H 2 (top right), C + (bottom left) and C (bot- 
tom right) abundances for the los models. Contours mark to- 
tal gas densities 20, 100, 500, 1000 and 2000 cm" 3 . C + and C 
abundance maps are clipped below 2 10" 4 cm -3 and 10" 4 cm" 3 , 
respectively. Lines of sight 0, 1, 2 and 3 on the ft(H 2 ) map refer 
to the discussion in the text. 
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Fig. 8. Same as Fig. but for CO (top left), CS (top right), 
CH (bottom left) and CN (bottom right). Abundance maps are 
clipped below 2 10" 5 cm" 3 for CO, 2 10" 10 cm" 3 for CS, 
2 lO" 8 cm" 3 for CH and 2 10" 9 cm" 3 for CN. 




n H [cm 3 ] 

Fig. 9. Abundances of H 2 , C + , C, CO, CH, CS and CN versus 
total gas density hr- Data points are averaged in the same hh 
bins as on Fig. [6] The vertical line marks the position of the 
average molecular transition where (fn 2 ) = 2 (ft(H 2 )/ft#) = 1/2. 



et al. 



(2008) are shown as dashed lines, and on the CO panel, 
we also plot actual data points from that same paper. 

Regarding CO, our data shows a significant deficit around 
A/"(H 2 ) ~ 10 20 cm -2 , which may be a general issue with PDR 
chemistry computations ([Sonnentrucke r et al~l |2007). However, 
the agreement with the observational fit gets better, both in val- 
ues and in slope, for ./V(H 2 ) > 2 10 20 cm -2 . The slope found 
is d[log#(CO)]/d[log#(H 2 )] - 5.2, while the observational fit 
yields 3.07 ± 0.73. On the longer lines of sight, we can see a 
"loop" structure which needs to be understood. It is easily iden- 
tified with lines of sight between Y = -24.5 pc and Y = -23 pc 
(positions to 3 on the H 2 and CO maps from Figs. [7] and [8}. To 
be more accurate, following the loop clockwise corresponds to 
scanning lines of sight parallel to X from to 3, with turnovers 
at positions 1 and 2. Considering the H 2 and CO maps alongside 
Fig. 11 helps to understand this "loop". Firstly, lines of sight 
from to 1 basically intercept just one dense structure with both 
H 2 and CO, so that we have a very similar behaviour to that of 
the short lines of sight parallel to Y. Secondly, lines of sight from 
1 to 2 pass through many less dense CO structures, with a lot 
of molecular hydrogen in between, leading to the sharp drop in 
the N(CO) vs. N(H 2 ) relation. At a given N(CO), the difference 
AAf(H 2 ) between H 2 column densities in the branches 0-1 and 
1-2 thus represents the "dark gas" (see subsection |5.4| ). Finally, 
lines of sight from 2 to 3 barely have any CO and contain less 
and less H 2 as we approach position 3, where we return to a sit- 
uation similar to that at position 0. The split between the two 
branches is definitely related to the mean gas density on the line 
of sight, the upper branch having (n H ) ~ 250 cm -3 , the lower 
one having (hh) ~ 150 cm -3 . The overall deficit in CO suggests 
that we may not shield it enough from the ambient UV field, a 
hypothesis which we discuss in section [6] 



The behaviour of Af(CH) with respect to N(U 2 ) (Fig. 11 



middle pan el) is in remar kable agreement with the observational 
data fit by 



Sheffer et al. 



(2008). Only for N(U 2 ) < 10 2U cm 



there a discrepancy, and it is irrelevant, as there are no detections 
in that range, only upper limits. 

bottom panel), we typically have a 



Concerning CN (Fig. 1 1 



factor 10 deficit when comparing to observations. It appears that 
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5.3. Comparison of the los and uniform models 



n(CH)/n(H 2 
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To assess the effects of taking into account density fluctuations, 
as opposed to the assumption of a homogeneous medium, usu- 
ally made when modelling PDRs, we compute the column den- 
sities of CO and H2 derived from the uniform models, and plot 
them on Fig. [12] The behaviour at low A/XH2) is very similar 



to that in the los models (Fig. [TT]), and the data suffer from 
the same deficit in CO around A/XH2) ~ 10 20 cm -2 . However, 
there is a definite difference at higher H2 column densities, as 
the slope of the relation between both column densities is much 
steeper, d[log #(CO)]/d[log N(U 2 )] - 14, than in the los mod- 
els and in observational data fits. This break occurs later, at 
Af(H 2 ) ~ 5 10 20 cm -2 , and applies to the short lines of sight 
(parallel to Y and parallel to X between positions and 1) where 
there is essentially one structure in both H2 and CO. However, 
the maximum column densities reached are slightly less than in 
the los models by a factor ~ 3, for both CO and H2. This shows 
the importance of taking into account density fluctuations along 
the line of sight when modelling PDRs. 

For CH, the behaviour is very similar in the uniform and 
los models, although here also the maximum column densities 
reached are slightly less in the uniform models, by a factor ~ 2. 
The observational fit is recovered at somewhat higher H 2 column 
densities (2 10 20 cm -2 instead of 10 20 cm -2 ), and the scatter of 
data points is a bit larger. 

Finally, regarding CN, if we ignore the deficit already seen in 
the los models, we reach the same conclusions : In the uniform 
models, a break in the slope occurs at higher H2 column densities 
N(U 2 ) ~ 6 10 20 cm" 2 (instead of N(U 2 ) ~ 4 10 20 cm" 2 in the 
los models), the slope is definitely steeper in that high-column- 
density regime, and the maximum column densities reached are 
smaller, here by a factor ~ 4. 



n(CS)/n(C) 




0.00090 
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0.00045 
0.00030 
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Fig. 10. Abundance ratios n(CH)/n(H 2 ) (top), n(CN)MCO) 
(middle) and n(CS)/n(C) (bottom) in the vicinity of the total gas 
density peak, located at (X ^ -4.7 pc, Y ^ -24 pc). Contours 
mark total gas densities n H of 100, 200, 300, 400, 500, 1000 and 
2000 cm" 3 . 



the reaction rate coefficient for the CN + N — > C + N2 reaction in 
the KIDA databas^jmay be too high (E. Roueff, private comm.), 
which could partly explain that deficit. 



5.4. Dark molecular gas fraction 

To estimate the amount of molecular gas not seen in CO - the so- 
called "dark gas" (Gren ieret al. 1 12005[ |Planck Coll aboration 
201 lj) - in our simulation, we need to compute the CO line emis 



sion and compare its spatial distribution with that of H2. To that 
effect, we focus on the cloudlet at (X ^ -4.7 pc, Y - -24 pc), 
where the gas density peak is found, and assume a very crude 
cylindrical cloud model by replicating the density maps along 
the Z axis over a line-of- sight L = 1 pc, which roughly corre- 
sponds to the cloudlet's extent in the (X, Y) plane. This effec- 
tively yields column density maps which we can use with the 
RADEX radiative transfer code ( [Van der Tak et al. ||2007| ) to ob- 
tain emission maps in the CO (/ = 1 — > 0) rotational transition 
line at 115.271 GHz and in the [CII] fine structure transition line 
at 158 yum. To be precise, at each position (X, 7), we treat radia- 
tive transfer along Z in a plane-parallel slab geometry. RADEX 
works with the escape probability formalism ( Sobolev 



1960), 



which requires specification of the line width A V. We estimate it 
to be cr 3D / V3 = 1 km.s -1 , where 0-3D = 1.8 km.s -1 is the total 
gas velocity dispersion listed in Table [T] Indeed, for any species 
a with molecular weight [i a (j^co = 28 and yUc+ = 12), the ra- 
tio of thermal to one-dimensional turbulent velocity dispersions 
reads 

&Ua) _ / T 



2.1 x 



'3D 



(270 k) 



http : //kida . obs . u-bordeauxl . fr/ 



In the cloudlet under study, T < 100 K, so the above ratio is typ- 
ically < 0.03 for CO and < 0.06 for C + . It is thus reasonable to 
take AV = cr 3D / V3 for all RADEX runs. The code also requires 
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specification of the gas kinetic temperature and the densities of 
collisional partners (H 2 for CO; H 2 , H and electrons for [CII]), 
which we get from the PDR code outputs. RADEX is thus run 
on every line of sight parallel to the Z axis, and results are com- 
bined into a CO (/ = 1 —> 0) emission map and a [CII] emission 
map, both shown on Fig. 13 

Between the solid and dotted contours is the "dark molecu- 
lar gas" region where hydrogen is predominantly in its molec- 
ular form but CO emission fails to detect it. We assume a 
detection threshold Wco = 0-4 K.km.s -1 consistent with the 
noise level in e.g. the CO survey of Taurus by Goldsmith et 



al. H2008) . On the other hand, that same gas can definitely be 



traced in the [CII] line, as it has a typical integrated emission of 
0. 4-0.8 K.km.s -1 , while the sensitivity quoted by Velusamy et 
al~| ( |2010t for the GOTC+ key program is ~ 0.1 - 0.2 K.km.s -1 . 

The "dark molecular gas" fraction associated with this 
cloudlet can be estimated by taking one-dimensional cuts par- 
allel to the Y axis going through the CO emission region. Along 
such a cut, which is parametrized by X, we define YnQC) and 
Y\ (X) as the boundaries of the computational domain^] (sharp 
transition from grey to white on the panels of Fig. 13 ), and we 
note Wco(X) the region where the integrated emission of CO 
(J = 1 — > 0) exceeds the detection threshold Wqo (region en- 
closed by the solid contour on Fig. 
gas" fraction as 



13 ). We then define the "dark 



/dg® = 1 " 



J w 



(H 2 )d7 

W C0 (X) 
Yi(X) 

(H 2 )d7 

Yo(X) 



Fig. 14 shows this fraction as a function of the position X of 
the one-dimensional cut. Obviously, /dg = 1 when the cut does 
not pass through the CO emission region, and /dg < 1 when 
some of the H2 is traced by CO. We find that in this cloudlet, 
at least 20% of H2 is not traced by CO, even at the peak of the 
gas density. To get a mean fraction of dark gas in this cloudlet, 
we average /dg over the range of X coordinates where CO is 
seen (i.e. /dgPO < 1), weighted by the total gas column density 
N H . This y ields fpc = 0-32, which is somewhat higher than the 
findings of Velusamy et al. ( 2010| ), who identified 53 "transi- 
tion clouds" with both Hi and iZ CO emission but no 13 CO, and 
found that ~ 25% of H 2 in these clouds belong to an H 2 /C + layer 
not see n in CO. However, the scatter in ob servational values is 
large ( [Grenier et al~| |2005[ |Abdo et al] |2010| ), so the small 



discrep ancy is no cause for al arm. Another possible comparison 
is with Wolfire et al. ( 2010| ), who constructed spherical mod- 
els of molecular clouds to study the dark gas fraction, which 
they define in a similar way, except that the boundary of their 
CO region is specified by the condition of unit optical thick- 
ness at the line center, tqo = 1. In our clump, that isocontour 
is very close to our own condition 7 C o = Wco> as can be seen 
on Fig 13 Computing the average dark gas fraction with their 
condition yields /dg = 0.36, which is quite close to their re- 
sults /dg ~ 0-25 - 0.33. There is a notable difference between 
their models and ours, however, since our cloudlet has a mass 
~ 9.5 M Q inside the CO region, while |Wolfire et al. | ( |2010| ) study 
GMCs with masses in the range 10 5 M Q to 3 10° M Q . Their im- 
pinging UV field is also notably higher (x = 3 - 30). 



6. Discussion and summary 

6.1. Illumination effects 

In this paper, we bypass the one-dimensionality of the PDR code 
by combining runs in two orthogonal directions, taking, at each 
grid point, the chemical composition corresponding to maxi- 
mum radiation energy density E. It should be noted that this 
choice is questionable : if a position is shielded from radiation 
in almost every direction but for one small hole, the illumina- 
tion at this location resulting from our procedure is too high. 
As this means forming less molecules, we wish to estimate if it 
might account for some of the CO deficit seen on Fig. [TT] To 



do so in a simple way, we compute the chemical composition in 
the opposite assumption, i.e. based on the criterion of minimum 
local illumination. The result is plotted on Fig 15 and shows 
how indeed this helps recovering observed CO column densi- 
ties for Af(H 2 ) > 10 20 cm -2 , with a consistent scatter. Below 
Af(H 2 ) - 10 20 cm -2 , a significant CO deficit remains, however. 

Obviously, this choice of minimum local illumination is also 
unphysical, and the reality must lie somewhere in between. A 
physically better, but more computationally intensive method is 
being pursued and will be presented in a future paper. 



6.2. FGK approximation 



Our study uses the |Federman et al. (1979]) (FGK) approximation 
to compute self- shielding. This may underestimate the shielding 
of CO by molecular hydrogen lines, so we perform a few runs 
of the PDR code using exact radiative transfer ( [Goicoechea & 
Le Bourlot 1 12007 ). We do this on some lines of sight for which 



N(H 2 ) ~ 10 /u cm -2 , to see whether this helps fill the CO deficit 
in that region. It turns out that the CO column densities so ob- 
tained are indeed higher than those found in the FGK approxi- 
mation, but by a factor < 2, which is not enough to explain our 
CO deficit. As the computational time is on the other hand in- 
creased by a factor ~ 5 - 6, we feel that this approach is not to 
be pursued. 



6.3. Steady-state assumption 

In this study, the simulation cube is taken as a static background, 
under the assumption that timescales for chemistry and photo- 
processes are much shorter than those of the MHD simulation. 



From the analysis performed by |Le Petit et al. | ( [2006| ) on uni- 



3 The dependence of / DG on the boundaries of the computational do- 
main is necessarily small, as there is little mass at low densities. 



form density PDRs, it appears that timescales for H 2 photodis- 
sociation at the edges of a cloud are ~ 1000/^ yr, where x is me 
FUV radiation strength in units of the |Draine~| ([T978) field. In the 
analysis presented, we choose x — 1 so mat me corresponding 
timescale is about 1000 yr. 

Estimating timescales for the MHD simulation is more dif- 
ficult, because what we're actually interested in is the time over 
which structures remain coherent, and we do not have access 
to this information due to the Eulerian nature of the simulation, 
which makes it impossible to confidently identify structures. For 
a rough estimate, we may consider the overall crossing time 
Across = L/Vx - 2.4 Myr, but this does not correspond to the 
time over which gas is mixed by turbulence at a given scale. For 
this, we may use the velocity dispersion cr 3D listed in Table [T] 
within the observational clump, whose size is about 2 pc. This 
yields a dynamical timescal e Td yn -1.1 M yr, which is very sim- 
ilar to the values quoted by Wolfi re et al. | ( |2010| ) to validate the 
steady-state assumption in their models. 
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The chemical timescale is that of the formation of molecular 

( [2007] ) using nu- 



hydrogen. As shown by |Glover & Mac Low 



merical simulations of decaying ISM turbulence that include a 
simplified chemical network, the formation timescale for H2 in 
turbulent magnetized molecular clouds is T C h em -1-2 Myr. It is 
therefore of the order of the estimated dynamical timescales in 
our simulation, so that our steady-state assumption seems only 
marginally valid. However, H2 is formed in dense regions and 
transported in the entire volume through turbulent motions, so 
we may be safe assuming steady state, provided we consider a 
late enough snapshot. Indeed, if H2 starts forming when the con- 
verging WNM flows collide near the midplane (r co n - 1 Myr), 
and if it is fully formed and transported in the entire volume after 
Tchem + Tdyn + T cross , this requires taking a snapshot timed at no 
earlier than 5.4 Myr, which is the case here it = 7.35 Myr). We 
conclude that our steady-state assumption is a legitimate one. 

6.4. Warm chemistry 

It should be noted that chemistry is here driven by UV radiation 
only, but that there is an important pathway for the formation of 
many molecular species, which is warm chemistry in turbulence 
dissipation regions (TDR), studied by Joulai net al. | ([1998) and 
|Godard et al. ( 2009). In particular , the CO abundances found in 



the models by Goda rd et al. ( |2009| ) are larger than in correspond- 



ing PDR models, sometimes by almost an order of magnitude. 
More generally, [Godard et al. | ( |2009 ) argue that observed chemi- 
cal abundances are on the whole well reproduced if dissipation is 
due to ion-neutral friction in sheared structures -100 AU thick. 
A TDR post-processing of our MHD simulation is in the works, 
to compare both types of chemistry. 

6.5. Summary 

We have presented a first analysis of UV-driven chemistry in 
a simulation of the diffuse ISM, by post-processing it with the 
Meudon PDR code. Our results show that assuming a uniform 
density medium when modelling PDRs leads to significant er- 
rors : in the case of CO, for instance, the maximum column 
densities found with this simplistic assumption are a factor ~ 3 
lower than those found using actual density fluctuations. The 
slope of the H 2 -CO correlation at N(H 2 ) > 5 10 20 cm" 2 is also 
a factor ~ 3 higher than in the more realistic case, and there- 
fore much less in agreement with observations. A second re- 
sult of our study is that, in the densest parts of the simulation 
(n H > 10 3 cm" 3 ), some 35% to 40% of the molecular gas is 
"dark", in the sense that it it not traced by the CO(/ = 1 — > 0) 
line, given current sensitivities. It is however detectable via the 
[CII] fine structure transition line at 158 yum. As a side result, 
we find that the simplified cooling used in the MHD simulation 
by [Hennebell e et al. | ( |2008| ) yields gas temperatures in reason- 
able agreement with those found using the more detailed pro- 
cesses included in the PDR code. 
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Table A.l. Parameters used in the PDR code, as input in the . in 
files. 



Parameter Description 



Appendix A: Using the Meudon PDR code with 
density profiles 

This appendix is meant as an introduction to using the Meudon 
PDR cod^jwith fluctuating density profiles. For more detailed 
presentations of the code, the rea der is referred to|Le Bourlot et 
|ai~| ( [T999l ); |Le Petit et al. [ ( [2006] ); [Gonzalez Garcia et al. | ( |20Q8| ). 

The code requires two input files : a .pfl file listing visual 
extinction Ay, temperature T (in K) and total gas density nn (in 
cm -3 ) along the line of sight, and a . in file supplying th e pa- 
rameters of the run to perform. These are listed in Table |A.1| 
and some of them require a short comment : 

- modele is the basename chosen for output files. For the los 
models, we use a name of the generic form los.zZ_yY_xX m -X p 
or los-zZ_xX_yY m -Y p , reflecting the position of the extracted 
profile in the cube. Note that Z is a constant throughout this 
paper, as is evident from Fig. [2] For the uniform models, 
we use names of the generic form uniform_zZ_y7_xX m -X /? or 
uni f orm_zZ_xX_y Y m - Y p . 

- ifafm is the number of global iterations to use. As Le Petit et 



al. ( 2006) point out, for diffuse clouds (Ay < 0.5) proper conver- 



gence may require up to 20 iterations, so we select if afm=2® for 
all of our models. For information, 415 of the 447 profiles have 
total A y < 0.5. 

- Avmax is that same total visual extinction through the cloud, 
which is simply the last Ay value in the . pfl file. 

- We set the density densh, temperature tgaz, and pressure 
presse=denshxtgaz parameters to the average value^] for 
each profile. They are not used by the code with a density- 
temperature profile, but they are used, in the reference uniform 
models, as initial guesses for thermal balance computation. 

- radm and radp specify the strengths^ and^ of the incident 
radiation field in units of the ISRF, respectively on the left and 
right sides of the profile. For the runs described in this paper, we 
use* m =x P = 1. 

- fprofil specifies the .pfl density-temperature profile file. 
For consistency, we use the same naming scheme as for modele. 

- vturb is the "turbulent velocity dispersion". It does not include 
thermal dispersion, so we take it to be the standard deviation of 
the line-of-sight velocity within each structure, noted crx and cry 
in Tables |2] and |3j respectively. 

-ifisob is a flag specifying whether to use a density profile. 
For the los models, we therefore set ifisob=l, to enforce the 
use of a density-temperature .pfl file. However, since thermal 
balance is solved (ieqth=l), the temperature values in the file 
are only used as initial guesses. For the uniform models, we set 
if isob=® to use a constant density (specified by densh). 



Appendix B: Post-processing of raw outputs 

B.1. Resampling 

Outputs of the PDR code are FITS data files and XML descrip- 
tion files, and we use dedicated scripts to extract specific quan- 
tities into plain text files for subsequent analysis. Among the 
quantities retrieved are the distance d from the surface of the 



modele 

chimie 

ifafm 

Avmax 

densh 

F_ISRF 

radm 

radp 

srcpp 

d_sour 

fmrc 

ieqth 

tgaz 

ifisob 

fprofil 

presse 

vturb 

itrfer 

j f gkh2 

ichh2 

los_ext 

rrr 

cdunit 
alb 

gg 

gratio 

rhogr 

alpgr 

rgrmin 

rgrmax 

FJDUSTEM 

iforh2 

istic 



Basename for the output files 

Chemistry file 

Number of global iterations 

Integration limit in A v 

Initial density (cm -3 ) 

ISRF expression flag 

ISRF scaling factory 

ISRF scaling factory 

Additional radiation field source 

Star distance (pc) 

Cosmic rays ionisation rate (10 -17 s -1 ) 
Thermal balance computation flag 
Initial temperature (K) 
State equation flag 
Density-Temperature profile file 
Initial pressure (K.cm -3 ) 
Turbulent velocity (km.s -1 ) 
UV transfer method flag 
Minimum J level for FGK approximation 
H + H 2 collision rate model flag 
Line of sight extinction curve 
Reddening coefficient R v = A V /E B _ V 
Gas-to-dust ratio C D = N H /E B _ V (cm -2 ) 
Dust albedo 

Diffusion anisotropy factor (cos 6) 
Mass ratio of grains / gas 
Grains mass density (g.cm -3 ) 
Grains distribution index 
Grains minimum radius (cm) 
Grains maximum radius (cm) 
DUSTEM activation flag 
H 2 formation on grains model flag 
H sticking on grain model flag 



Value 

see appendix^ 
chimie08 a 
20 

see appendix A 
see appendix ~K 
l b 

see appendix A 
see appendix ~K 
none . txt 
C 
5 

l d 

see appendix A 
see appendix ~K 
see appendix ~K 
see appendix A 
see appendix A 
e 


2 f 

Galaxy g 
3.1 h 

5.8 X 10 21 g 
0.42 « 
0.6 « 
0.01 § 
2.59 1 
3.5§ 

3 x 10 -7 8 

3 x 10 -5 g 

j 

k 

4 1 



This chemistry file does no t include deuterated species 



Mathis et al. 



( |1983| > and |Black |(l994| 



Uses expr essions based on 
rather than |Draine ] ( |1978| ). 
This means that no additional star is present. 
Thermal balance is computed for each point in the cloud, using the 
tempera tures given by the . pfl fi le as initial guesses. 
Use the Feder man et al. | {1979) (FGK) approximation for the H 2 
lines in the UV. 

Use the values compiled b y Le Bourlot et al. ( |1999| > with reactive 



collisions from Schofield (1967). 

See Table 4 of |Le Petit et al. | i |2006| ), which quotes val ues from 
Fitzpatrick & Massa (1990 ) for the extinction cur ve, |Bohlm~"| 
( |1978| ); IRachford et al. |p002J> for C D ,IMathis | ( [T996]) for the dust 
albedo, Weinga rtner & Draine | ( |2001| ) for (cos iff), |Mathis et al.~| 
(1977) for the grain size distribution parameters. 
See for inst ance |Cardelli et al.l(|l989|> 

Taken from Gonzalez Garcia et al. ( 20081. 



We use version 1.4.1 of the PDR code, with a fixed H 2 formation 
rate R f = 3 10 -17 V^/100 K. 
5 i.e. density-weighted averages for temperature and pressure. 



Do not couple to the DUSTEM code (Compieg ne et al. | 12010). 
Energy released by H 2 formation on dust grains is equally split be- 
tween grain ex citation, H 2 kinetic e nergy and internal energy. See 
section 6.1.2 in Le Petit et al. ( 2006 ) for details. 
See Appendix E5 in |Le Petit et al. | < [2QQ6] >. 



structure, visual extinction Ay, proton column density Nh, tem- 
perature r PDR , proton density n H , pressure p, ionization fraction 
x e , and abundances n(a) of 99 chemical species. 

As the PDR code does its own mesh refinement to better 
solve for the H/H 2 transition, these quantities are sampled irreg- 
ularly. Consequently, we resample outputs on the same regular 
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grid as the MHD simulation, using a simple linear interpolation 
method. This allows us to build raw maps for all output quanti- 
ties from PDR code runs along the X and Y directions. 

B.2. Missing data 

For reasons that are unclear, a few^] of the 894 PDR runs do 
not complete successfully on the EGEE grid. To supplement the 
missing data, we interpolate along the perpendicular direction. 
Consider the ensemble of runs along the X direction : completed 
runs yield quantities F X (X, Y), and if a run is missing at coor- 
dinate Y = Yq, we supply F X (X, Yq) by linearly interpolating 
G(Y) - F x (Xq, Y) at constant X = Xq. This yields a satisfactory 
completion of the raw data. 

B.3. Combination ofX and Y runs 

We then proceed to the combination of d ata fr om runs along the 
X and Y directions, as described in[3] Fig. A.l shows the "illumi- 
nation mask" computed by comparing the local radiation energy 
densities E x and E Y output by the PDR code in los models par- 
allel to the X and Y directions, respectively. This mask is then 
used to build a single data array for each quantity F output by 
the PDR code, at each grid point (X, F), according to the rule : 



(n H ) [cm 3 ] 



F x 
F Y 



if 
if 



E X > Ey 

E x < E Y 



This helps to reduce the shadowing artifacts due to the one- 
dimensionality of the PDR code, while ensuring element con- 
servation in each grid cell, and yields the final maps that are 
analyzed and discussed in the main body of the paper. In the 
discussion (section [6]), we also make use of the inverse choice : 



F = 



F Y 
F x 



if 
if 



E x > E Y 
E x < Ey 




Fig. 11. Column densities of CO (top), CH (middle) and CN 
(bottom) versus column densities of H2, in the los models. 
Circles correspond to lines of sight parallel to Y and squares to 
lines of sight parallel to X. Their colours reflect the mean gas 
density (nn) on these lines of sight. Plus signs on the top (CO) 



plot stand for observational data points (Sheffer et al. 2008 
Crenny & Federman 1 12004[ [Pan et al. | 2005[ |Lacour et al. 
2005j I Rachford et al. \ |2002l |20091|Snow et al. 1 12008|) . The 



dashed lines are power-law fits from Sh effer et al. | (|2008 ). The 
lines of sight parallel to X marked 0, 1, 2 and 3 on the top panel 
are the same as on Figs. [7] and [8] 



6 Namely, for the los models, 16 out of 291 along Y and 3 out of 156 
along X; for the uniform models, 8 out of 291 along Y and 6 out of 156 
along X. 
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Fig. 13. Synthetic emission maps in CO(/ = 1 — > 0) (left) and 
[CII] at 158 yum (right). The solid contour marks the assumed 
0.4 K.km.s -1 detection threshold for CO, the dashed contour 
marks the position of line center optical depth tqo = 1, and 
the dotted contour marks the position of the molecular transition 
[ Cm -3] fn 2 = 1/2. Grey areas are outside of the computational domain. 
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Fig. 14. Fraction of "dark gas" (solid line) along one dimensional 
cuts parallel to the Y axis going through the cloudlet at (X 
-4.7 pc, y -24 pc). Also shown are the fraction of "dark gas" 
computed using the definition by Wol fire et al. | ( |2010| ) (dashed 
line), and a normalized profile of the total gas column densities 
Nh along the same cuts (dash-dotted line). 
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Fig. 12. Same as Fig.[lj]but for the uniform models. 
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Fig. 15. CO column densities versus H2 column densities, in the 
los models when combining data based on a criterion of mini- 
mum local illumination. Symbols are the same as on Fig. [TT] (top 
panel). 
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Fig.A.l. Illumination mask computed by comparing at each 
point (X, Y) the local radiative energy densities E x and E Y out- 
put by the PDR code along the X and Y directions, respectively. 
Regions where E Y > E x are marked in black and regions where 
E x > E Y are marked in white. Contour lines of equal total gas 
density are overlaid at 20, 100, 500, 1000 and 2000 cm" 3 . Grey 
areas are outside of the computational domain. 
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